Ray marching is an alternative to analytic methods for discovering intersections between a light transport ray and the scene. The structure of the ray marching algorithm makes it amenable to rendering scenes constructed from operations on solid objects, amorphous surfaces, and participating media, which have no explicit surface.

Analytic intersection is extremely efficient. Renderers based on it can render simple scenes interactively on a CPU, which is great for learning and experimentation. Ray marching is less efficient. So, this chapter also introduces Graphics Processing Unit (GPU) programming to enable you to build a real-time ray marching renderer.

© 2014 Morgan McGuire
A reference scene defined by implicit surfaces ray traced at 60fps on
GeForce 650M using the techniques described in this document.

© 2013 Pablo Roman Andrioli
An impressive open source demo of real-time ray-traced implicit surface fractals
by Kali (Pablo Roman Andrioli), implemented in about 300 lines of GLSL.

1. Intersection AlgorithmsGPU implicit signed distance

Most computational graphics rendering methods are based on ray-surface intersection. The three major ray-surface algorithms are:

• Analytic ray intersection
• Numerical intersection by ray marching
• Rasterization with a depth buffer

Analytic ray intersection and ray marching use implicit surface definitions. For example, the unit-meter sphere is implicitly defined by the set of points satisfying

That is, the zero-crossing of the function h∈R3m→Rm measuring signed distance from the surface of the sphere.

Analytic ray-surface intersection gives a closed-form solution by solving for the distance parameter t for which the explicit equation of a ray,

satisfies the implicit equation of the scene, in this case, a sphere. While analytic solutions exist for some cases, such as ray-sphere and ray-triangle intersection, there is no general analytic solution possible for arbitrary ray-surface intersection.

Ray marching is an algorithm for numerical estimation of arbitrary ray-surface intersections. It works for arbitrary implicit surfaces (although most of the optimizations that make it practical for real-time use require a surface definition with certain properties).

The simplest form of ray marching begins at the ray origin and samples the implicit surface at regular intervals, terminating when the function crosses zero. This is of course both slow and innaccurate. Fortunately, any numerical root-finding strategy, such as the secant method or Newton-Raphson iteration, can be employed to accelerate it or improve its accuracy. For example, choosing the step size between successive samples based on the derivative of the surface gives faster convergence and more accuracy. Binary search near the zero crossing is a good way to identify the intersection with arbitrary accuracy for continuous implicit surface definitions.

Rasterization with a depth buffer is an amortization of analytic ray intersection over many rays that relies on both explicit and implicit surface definitions. It processes many rays for each surface, instead of many surfaces for each ray. Because of this, it relies on an auxillary image-space structure called the depth buffer to keep track of the currently-known closest surface location to the camera at each image sample. (There are also rasterization methods that perform sophisticated clipping of primitives against each other to avoid this structure, however their increased run time cost puts them out of favor today.)

Rasterization is extremely efficient for triangle meshes on massively parallel hardware, and algorithms are known for some other primitives such as lines, rectangles, and disks. However, because it requires both analytic implicit and explicit equations, rasterization is highly constrained to such primitives.

2. Marching

Consider an example of solving for the root of the equation that gives the intersection of a ray P+ω^t and a general sphere with center point C and radius r,

The simplest implementation of ray marching iterates t through the domain of function g in fixed Δt steps:

1. 1	def findRootNaive(g∈Rm→Rm)∈Rm:

2. 2	Δt∈Rm=0.01m

3. 3

   t∈Rm=0m

4. 4	while g(t)≠0m:

5. 5

   t=t+Δt

6. 6

   return t

This “marches” through every point X(t)=P+ω^t on the ray until it finds a t that yields an X on the surface of the sphere.

Of course, this naive ray marching scheme is unlikely to exactly land on the exact value of t for which g(t)=0m, so the code will almost always loop forever. A better solution enforces a maximum distance bound and seeks the first value at which g crosses from positive to nonpositive, at which point the ray must have passed through the surface:

1. 1	def findRoot(g∈Rm→Rm, maxDistance∈Rm)∈Rm:

2. 2	Δt∈Rm=0.01m

3. 3

   t∈Rm=0m

4. 4	while (t < maxDistance) and (g(t)>0m):

5. 5

   t=t+Δt

6. 6

   return t

We can further refine this approach using binary search within the final fixed interval [tn−1,tn] in which g crossed from positive to nonpositive values.

Ray marching with fixed steps has three advantages over analytic ray intersections:

• Very simple to implement
• Works for any implicitly defined surface
• Naturally extends to integration of scattering in participating media, such as fog

Ray marching with fixed steps has two obvious disadvantages over analytic ray tracing. The first is that it can miss the intersection for a non-convex shape. The second is that it takes at least linear time in the distance to the intersection for rays that do intersect the surface and linear time in maxDistance for rays that do not intersect. Analytic ray tracing is guaranteed to not only find the intersection, but do it in constant time.

We can mitigate both the performance disadvantage and the missed intersection disadvantage of fixed-step ray marching with a simple improvement that only slightly constrains the surface description.

3. Distance Estimators

We originally defined g(t)=h(P+ω^t)=0m as the surface. The ray marching algorithm implicitly required only that g(t)>0m at the ray origin P and g was continuous. Now, let us require the surface to be an embedded manifold without boundary, i.e., the surface of a solid object, and refine our constraints so that g and h are signed distance estimators: require that g(t)=h(P+ω^t) is less than or equal to than the true distance to the surface from point X(t)=P+ω^t.

Because we're working with surfaces in three dimensions and will wish to cast rays from many different origins (e.g., primary rays, shadow rays, mirror rays), it will be convenient to define a distance estimator on three-space instead of the ray parameter.

We now formally define function h(X)=g(P+ω^t) be a signed distance estimator for a surface. For the sphere, one function that satisfies this is simply the original implicit equation: hsphere(X)=||X−C||−r.

The h(X)=0m surface is sometimes called the isosurface, although technically it is the h=0m isosurface. Fluids are often modeled with a large number of particles, which are collectively rendered by forming a single distance estimator that is the sum of a large number of so-called blobby distance estimators. In this case, the term level set is sometimes used to describe the surface.

4. Sphere Tracing

Given a signed distance estimator, h(X), we now have an upper bound on how far is safe to step to avoid marching the ray over the intersection. Because h returns a bound on the distance in any direction, we can move at least that far along the direction ω^ of the ray. In other words, if we consider a sphere of radius h(X) of possible locations to visit from point X, that sphere at most touches (i.e., is tangent to) the primitive (...which is also a sphere in our example) that we're tracing against, and may not even touch it.

That new virtual sphere of possible locations to visit can't possibly intersect the primitive more deeply because of the distance estimator guarantee. So, we can jump to testing the next point at the surface of the virtual sphere. This method for advancing the ray is called sphere tracing [Hart1996Sphere] . Note that the “sphere” in the name refers to the sphere within which the algorithm can safely march from its current location; the method can trace any primitive, and intersects the primitives with rays, not spheres.

1. 1	def sphereTrace(P∈R3m,w∈S2,h∈R3m→Rm,maxDistance∈Rm)∈Rm:

2. # How close we want to get to the surface

3. 2	ϵ∈Rm=0.0000001m

5. # Enforce this minimum step through space

6. 3	minStep  ∈Rm = 0.0001m

8. 4

   t∈Rm=0m

9. 5	Δt∈Rm=h(P+ω^t)

10. 6	while (t < maxDistance) and (Δt > ϵ):

11. 7	t=t+max(Δt,minStep)

12. 8	Δt=h(P+ω^t)

14. 9

    return t

Why is iteration necessary? Doesn't h give the distance to the surface immediately? If h gave the exact distance to the surface, then the while loop would be unnecessary. However, recall that h is a conservative estimate of the distance. If the surface is 1m from a point X, then h(X) might return 0.1m. So, we must iterate until h(X) converges to zero, or at least very nearly so.

The minStep constant is necessary to ensure that the ray continues to advance even when asymptotically approaching a surface. It is particularly important for achieving reasonable performance in the face of an overly conservative distance estimator or at locations where the ray just barely passes by the surface. One could additionally impose an explicit iteration limit independent of the minimum step size. Keinert et al. [2014] recommend tracking and returning t corresponding to the closest point ever found to the ray instead of the marching end point when terminating after a fixed number of iterations.

As before, we can extend this root finding method to perform binary search to refine the end hit point. Because h is a signed distance estimator, we can even do so fairly efficiently. Once the ray passes through the surface we can re-approach it accurately from the inside.

5. Some Distance Estimators

For these definitions, let max(v⃗ )=max(xv,yv,zv), |v⃗ |=(|xv|,|yv|,|zv|), and let max(a⃗ ,b⃗ )=(max(xa,xb),max(ya,yb),max(za,zb)). The code samples are in GLSL using some common utilities and HLSL compatibility definitions:

Copy

1. 1	#define Vector2 vec2

2. 2	#define Point3 vec3

3. 3	#define Vector3 vec3

4. 4	#define lerp mix

5. 5	float maxComponent(vec3 v) { return max(v.x, max(v.y, v.z)); }

6. 6	float saturate(float x) { return clamp(x, 0.0, 1.0); }

Sphere

The sphere about point C with radius r:

Plane

The plane with point C and normal n^:

Box

The box with center C and vector b⃗  from the center to the first-quadrant corner:

Rounded Box

The hollow rounded box with center C, edge half-lengths in vector b⃗ , and rounding radius r:

Torus

The torus about the x-axis with centroid C, minor radius r, and major radius R:

Wheel

The “wheel” about the x-axis with centroid C, minor radius r, and major radius R is the same as the torus above, but with a non-Euclidean definition of distance:

Cylinder

The cylinder with centroid C, radius r, and half-height e:

6. Computing Normals

The gradient of any function on 3-space can be approximated by numerical differentiation:

If the distance estimator is not too conservative near the surface, or at least is conservative by the same factor in all directions, we can use this to find the surface normal.

For a point X that is near but not on the surface, normalizing the gradient by dividing through by its length gives the surface normal vector to the nearby surface.

Placing X exactly on the surface can be problematic for this computation, however. For example, for a sufficiently curvy surface, a particular ϵ value chosen may place all samples of h at another location on the surface, for which h=0m.

Reasonable four- and six-sample approximations of the normal are thus:

Copy

1. 1	Vector3 fastNormal(Point3 X) {

2. 2	const float e = 1e-4;

3. 3	float d = sceneDistance(X);

4. 4	return normalize(Vector3(

5. 5	sceneDistance(X + Vector3(e, 0, 0)),

6. 6	sceneDistance(X + Vector3(0, e, 0)),

7. 7	sceneDistance(X + Vector3(0, 0, e))) - Vector3(d, d, d));

8. 8

   }

10. 9	Vector3 normal(Point3 X) {

11. 10	const float e = 1e-4;

12. 11	const Vector3 u = Vector3(e, 0, 0);

13. 12	const Vector3 v = Vector3(0, e, 0);

14. 13	const Vector3 w = Vector3(e, 0, e);

15. 14	return normalize(Vector3(

16. 15	sceneDistance(X + u) - sceneDistance(X - u),

17. 16	sceneDistance(X + v) - sceneDistance(X - v),

18. 17	sceneDistance(X + w) - sceneDistance(X - w)));

19. 18

    }

7. A Simple GLSL Ray Caster

This section describes the ray casting portion of a simple GLSL ray marching tracer that runs on the GPU using the G3D Innovation Engine 10. See http://codeheartjs.com/examples/raytrace/ for a comparable CPU Javascript version that runs in a web browser, and http://codeheartjs.com/examples/fastraytrace/ for an optimized Javascript version. See my gpu-tracing-tutorial.zip for a working example of an analytic ray tracer in C++ and GLSL against which to compare this ray marcher's structure. You can inject the ray marcher directly into that program structure as well.

C++ Shader Launch

The GLSL ray marcher below uses rasterization to submit two giant triangles that cover the entire viewport, under a 2D projection matrix. It then launches the ray tracing shader kernel at each pixel, passing it the 3D camera's orientation and projection matrix. The C++ syntax for setting up this call in OpenGL (from App::onGraphics3D) is:

Copy

1. 1	rd->push2D();

2. 2	Args args;

4. 3	args.setUniform("cameraToWorldMatrix", activeCamera()->frame());

6. 4	args.setUniform("tanHalfFieldOfViewY",

7. 5	float(tan(activeCamera()->projection().fieldOfViewAngle() / 2.0f)));

9. // Projection matrix, for writing to the depth buffer.

10. 6	Matrix4 projectionMatrix;

11. 7	activeCamera()->getProjectUnitMatrix(rd->viewport(), projectionMatrix);

12. 8	args.setUniform("projectionMatrix22", projectionMatrix[2][2]);

13. 9	args.setUniform("projectionMatrix23", projectionMatrix[2][3]);

15. 10	// A cube map Texture

16. 11	m_skybox->setShaderArgs(args, "skybox.", Sampler::cubeMap());

18. 12	// Set the domain of the shader to the viewport rectangle

19. 13	args.setRect(rd->viewport());

21. 14	LAUNCH_SHADER("trace.pix", args);

22. 15	rd->pop2D();

GLSL Ray Setup

The GLSL code within the referenced trace.pix shader sets up the primary ray and then marches it through the scene. I've also rigged it to write a depth buffer value that is compatible with rasterization so that rasterized and ray traced primitives may be mixed in the frame buffer. Furthermore, this enables the use of standard post-processing tricks such as single-image depth of field, or even depth-buffer based ambient occlusion as a post process.

Copy

1. 1	#version 410 // -*- c++ -*-

3. 2	#include <g3dmath.glsl>

4. 3	#include <Texture/Texture.glsl>

6. // Input arguments from the C++ program

8. 4	uniform mat4x3 cameraToWorldMatrix;

9. 5	uniform float tanHalfFieldOfViewY;

10. 6	uniform float projectionMatrix22, projectionMatrix23;

11. 7	uniform TextureCube skybox;

13. // Output to the App::m_framebuffer (assume post-processed tonemap and gamma encoding)

14. 8	out Radiance3 pixelRadiance;

16. 9	void main() {

17. // Generate an eye ray in camera space, and then transform to world space

19. // Primary ray origin

20. 10	Point3 P = cameraToWorldMatrix[3];

22. // Primary ray direction

23. 11	Vector3 w = Matrix3(cameraToWorldMatrix) *

24. 12	normalize(Vector3((gl_FragCoord.xy - g3d_FragCoordExtent / 2.0) *

25. 13	Vector2(1, -1), g3d_FragCoordExtent.y /

26. 14	( -2.0 * tanHalfFieldOfViewY)));

28. 15	float distance = inf;

29. 16	pixelRadiance = traceRay(P, w, distance);

31. // Camera space z value

32. 17	float csZ = maxDist / w.z;

34. // Pack into standard OpenGL depth buffer format to make the result

35. // compatible with rasterization and post-processing.

36. 18	gl_FragDepth = (maxDist == inf) ? 1.0 :

37. 19	((projectionMatrix22 * csZ + projectionMatrix23) / -csZ);

38. 20

    }

GLSL Intersection

© 2013 Morgan McGuire
The result of the simple
GLSL ray marcher.

The actual tracing code is below. This example defines the scene to be a single, yellow sphere at the origin with a radius of 1m. To make it easier to orient the viewer when investigating the scene with the default camera, I also placed a cube map skybox in the background.

For a full ray tracer, we'd shade the intersection ray and possibly spawn additional rays for reflection, refraction, and shadows.

Copy

1. 1	float sceneDistance(Point3 X) {

2. 2	const Point3 C = Point3(0,0,0);

3. 3	const float r = 1.0;

4. 4	return length(X - C) - r;

5. 5

   }

7. // Trace the ray P + w * t.

8. // Returns true if something was hit.

9. // Sets L_o if the ray reached infinity or if it hit something before distance

10. // If something was hit, updates distance

11. 6	bool traceRay(Point3 P, Vector3 w, out L_o, inout float distance) {

12. 7	const float maxDistance = 1e10;

13. 8	const int maxIterations = 50;

14. 9	const float closeEnough = 1e-2;

15. 10	float t = 0;

16. 11	for (int i = 0; i < maxIterations; ++i) {

17. 12	float dt = sceneDistance(P + w * t);

18. 13

    t += dt;

19. 14	if (dt < closeEnough) {

20. 15	distance = t;

21. 16	L_o = Radiance3(1, 1, 0);

22. 17	return true;

23. 18	} else if (t > distance) {

24. // Too far; there is some known closer intersection

25. 19	return false;

26. 20

    }

27. 21

    }

29. // Reached infinity

30. 22	L_o = texture(skybox.sampler, w).rgb * skybox.readMultiplyFirst.rgb;

31. 23	return false;

32. 24

    }

This setup of using minimal OpenGL to launch a full-screen ray marcher is common in the modern demoscene. Many examples of artistically impressive demos and intros based around this technique can be found at pouet.net. Smaller single-shader examples with source code can be found at shadertoy.com and glslsandbox.com.

8. Operations on Distance Estimators

An advantage of modeling surfaces with distance estimators is that it is much easier to perform operations on whole shapes in that model than it is when they have explicit triangle representations or direct implicit definitions.

Operators on distance estimators are higher-order functions: they take functions as input and return functions as output. In a language such as Javascript or Scheme that has full support for closures, this can be implemented directly with functions. (As was done in the http://codeheartjs.com/examples/raytrace example from the previous section, which runs in a web browser.)

In a language such as Java or C++, classes and inheritance can mimic the design pattern. A simple substitution-compiled language like GLSL with provides none of these language features for abstracting computation. They can be emulated using macros and overloading or one can directly expand that code by hand, for example,

Copy

1. 1	float unionDistance(float d1, float d2) { return min(d1, d2); }

3. 2	float twoSphere(Point3 X) {

4. 3	return unionDistance(sphereDistance(X, Point3(-1, 0, 0), 1),

5. 4	sphereDistance(X, Point3( 1, 0, 0), 1));

6. 5

   }

9. Some Useful Operators

Union

Union selects the closer of two surfaces. This is equivalent to a set sum on the points (but not a distance sum of the estimators).

Intersection

Intersection selects the farther of two surfaces.

Subtraction

The set difference (subtraction) of two surfaces turns one inside out and then intersects its outside it with the other's inside.

Repetition

Any distance function can be tiled across space with period v⃗  across the dimensions.

Transformation

Transforming a shape is equivalent (from the reference frame of the distance estimator) to inversely transforming the points tested against it. In doing so, we need to be careful about how we have scaled space, however. To transform a shape by an invertible 4×4 rotation-translation-dilation (dilation means a uniform scale along all axes) matrix M,

Because rotation, translation, and non-zero dilation all have trivial inverses and determinants it is not usually necessary to apply the full inverse and determinant operations. For example scaling a primitive by scalar factor s>0 is simply:

Note that if a transformation of a simple primitive is desired, it may be more efficient to let C=(0,0,0) in the original primitive's definition so that it is centered at the origin and then perform a single transformation to the desired reference frame.

Blending

With some care to preserve the conservative property of distance estimation, shapes can be blended for a kind of smooth union. Hoffman and Hopcroft [1985] derive the general case of this, which is commonly referred to as a union using a “smooth min” (smin) function in place of a typical min function.

There are many smin functions, with a variety of performance and quality tradeoffs. Iñigo Quilez [2008] and I recommend this fast polynomial version:

Copy

1. 1	float smin(float a, float b, float blendRadius) {

2. 2	float c = saturate(0.5 + (b - a) * (0.5 / blendRadius));

3. 3	return lerp(b, a, c) - blendRadius * c * (1.0 - c);

4. 4

   }

10. Increasing Performance

Performance is always a concern in rendering. The performance of ray marching can be increased at the cost of image quality by increasing the minimum step size and accepted distance from the surface, and by decreasing the maximum iterations permitted and resolution. However, it can also be increased by algorithmic improvements that do not affect image quality. Instead, the cost is some of the elegant simplicity of the naive method. However, compared to analytic ray casting and rasterization, optimized ray marching remains remarkably straightforward and accessible to implement.

Over-Relaxation

Keinert et al. [2014] make an observation that allows decreasing the number of iterations before binary search by up to a factor of two. Hart's sphere tracing uses the distance estimator h to create a series of “unbounding” spheres as the ray approaches a surface. The h(X)=0m isosurface cannot lie within any of these spheres. Therefore, if the unbounding spheres at X and Y overlap, the surface cannot intersect the line segment XY.

This observation allows speculatively advancing the ray by Δt=k⋅h(X) for 1≤k≤2 under the following process. Let the new point be Y=X+ω^Δt. If h(Y)≥k⋅h(X), then there was no intersection on XY and the speculative advancement was conservative and can be accepted. Otherwise the ray can still be advanced from X by Δt=h(X). Keinert et al. call this over-relaxed sphere tracing (figure~\ref{fig:relaxation}. Naive sphere tracing advances the ray with only k=1, often discovering that the next point redundantly covers the interval just traversed.

© 2014 Keinert et al.
Top: sample points in white and ``unbounding spheres'' in purple from sphere
tracing. Bottom: Over-relaxed sphere tracing can speculatively take
larger steps (red) until it fails to find consecutive bounds
(yellow) and falls back to sphere tracing.

Over-relaxation is obviously most useful when the distance estimator in space is overly conservative for distance to the surface along the ray. There is a tradeoff between setting k close to 2 to gain the largest possible step and setting k close to 1 to reduce the number of times the speculative advancement fails. Because success doubles the step distance and failure doubles the number of evaluations of distance estimator h per iteration, over-relaxed sphere tracing can double or halve the performance compared to naive sphere tracing. Keinert et al. observe that k=1.2 gives a net 25% reduction in tracing time for the complex urban environment shown in the figure below, which appeared in the Mercury demoscene production The Timeless.

© 2014 Mercury
The Timeless demoscene production rendered with
sphere tracing optimizations to implicit surface ray marching

Bounding Spheres

When working with implicit surfaces, operations on distance estimators make it possible to build complex scenes by layering operations on simple primitives. Doing so offers great expressive power and programming elegance at a potentially high performance cost. This is because rays that pass nowhere near a surface must still evaluate the distance function for it.

The parse tree of a complicated set of operations on distance estimators inherently forms a scene graph tree. Performance can be increased by avoiding descending branches of that tree that do not affect the result. This kind of pruning is a common computer science operation and yields order-of-growth increases, unlike the kind of small constant factor that can be obtained by micro-optimizing within a single distance estimator.

The easiest way to prune the tree is to abort computation within a distance estimator when an early, inexpensive test determines that the ray is very far from the surface. A simple bounding sphere accomplishes this and correctly returns a still-conservative value. For example, let expensiveDE be a computationally-expensive distance estimator for an surface that is known to be entirely contained within the sphere with center C and radius r. A fast bounding estimator first tests how close the point X is to the bounding sphere. This must be a distance less than that to the true surface. It then only invokes the full, expensive estimator if X lies within the bounding sphere.

Copy

1. 1	float boundingDE(Point3 X, Point3 C, float r) {

2. 2	float t = length(X - C) - r;

3. 3	if (t > 0) { return t; } else { return expensiveDE(X); }

4. 4

   }

Square root is one of the most expensive operations to perform on a GPU. We can slightly increase the performance of the slow case of boundingDE without affecting the fast case by delaying the square root until the inside of the slow branch:

Copy

1. 1	float boundingDE(Point3 X, Point3 C, float rSquared) {

2. 2	Point3 v = X - C;

3. 3	float tSquared = dot(v, v) - rSquared;

4. 4	if (tSquared > 0) {

5. 5	return sqrt(tSquared);

6. 6

   } else {

7. 7	return expensiveDE(X);

8. 8

   }

9. 9

   }

It is also a good idea to branch to the full, expensive distance estimator if X is even close to the sphere. That is because the silhouettes of objects are the most expensive to trace using a numerical root finder. A ray passing by in a direction tangent to the surface appears to be very close to the surface to the omni-directional distance estimator, so the ray marching iterator takes very small steps, even though the ray is not actually approaching the surface quickly.

Introducing a false silhouette around each object through the bounding sphere can exacerbate this problem. We can separate the bounding radius from the test value to avoid that:

Copy

1. 1	float boundingDE(Point3 X, Point3 C, float rSquared) {

2. 2	Point3 v = X - C;

3. 3	float tSquared = dot(v, v) - rSquared;

4. 4	if (tSquared > rSquared * 0.2) {

5. 5	return sqrt(tSquared);

6. 6

   } else {

7. 7	return expensiveDE(X);

8. 8

   }

9. 9

   }

Reintroducing Analytic Roots

Geometric planes are are large surfaces that can cover a large amount of the screen. A ray caster following any ray near the plane will be constrained to march in steps no larger than the elevation of the ray above the plane, which may be quite small compared to the true distance to the next intersection with the scene.

We can obtain a good speedup for scenes that containh ground planes by performing the (efficient) analytic intersection with the plane and then using the ray marcher only for other surfaces. Spheres are a border case. The distance estimator is extremely fast and yields good estimates along the ray on the interior and far from the sphere. Yet, for rays that pass near the silhouette, an analytic solution is far faster than a ray marched one.

The figures below shows on top a rendering of a scene composed of implicit surfaces. On the bottom is a visualization of the number of invocations of the scene distance estimator at each pixel, where brighter pixels are more expensive. This scene uses the analytic ground plane and bounding sphere optimizations.

© 2015 Morgan McGuire
A reference scene of many different shapes defined by implicit surfaces.
The distance estimator and ray marcher was used for primary rays,
(soft) shadow rays, and ambient occlusion.

© 2015 Morgan McGuire
Visualization of the number of scene distance estimator invocations
at each pixel; black = 0, white = 20.

Note that in the above figures, silhouettes and surfaces seen at glancing angles are far more expensive than the interiors of objects, which are almost as efficient to trace as the ground plane itself. The screw-like shape and the heightfield are expensive because they have very conservative distance estimators. The silhouettes of shadows are also expensive--those are places where the ray for the visibility query to the light passes close to a surface and thus slows down as it approaches.

Other Optimization Strategies

Soft shadows, ambient occlusion, and in some cases antialiased pixels and depth of field can all be approximated efficiently using distance estimators because they track how close a point is to all surfaces in the scene. Under analytic ray tracing, these instead must be estimated using tens or hundreds of rays per pixel.

Once objects are surrounded in bounding spheres, traditional spatial subdivision data structures can be applied to them. For example, it is simple to compute a bounding volume hierarchy on the bounding spheres.

However, using a data structure on a GPU in graphics mode is a bit tricky because GLSL executes function invocation by inlining. Thus there is no true recursion for traversing a tree in a natural way. This also makes it hard to implement reflection and refraction at the same surface under Whitted ray tracing, since that creates a tree of rays. (Path tracing, photon tracing, and one of reflection or refraction without the other can all be implemented using loops.)

There are two ways to build an explicit stack to work around the lack of recursion. One is to build the stack global memory using OpenGL 4 image buffer operations. This has awkward syntax and can be slow. Another is to build the stack in local memory. This can be accomplished for recursive rays on recent hardware using code like the listing below. A similar explicit stack can be used to traverse data structures.

Copy

1. 1	struct RayStackFrame {

2. 2	Point3 origin;

3. 3	Vector3 direction;

4. 4	Color3 weight;

5. 5

   };

7. 6	RayStackFrame stack[MAX_STACK];

8. 7	int stackTop = 0;

10. 8

    ...

11. 9	stack[0] = RayStackFrame(X, w, Color3(1, 1, 1));

12. 10	while (stackTop >= 0) {

13. 11	Color3 weight = stack[stackTop].weight;

14. 12

    ...

15. 13	if (reflection) {

16. 14	// cast a recursive ray by pushing onto the stack

17. 15	// and attenuating by the magnitude of the BSDF impulse

18. 16

    ...

19. 17	++stackTop;

20. 18

    }

21. 19	L_o += weight * L_scattered;

22. 20

    }

Beware that local memory is significantly slower than register memory (where other local variables are stored), although significantly faster than global memory. Local memory is frequently implemented as a portion of the L1 cache. Also beware that on older hardware without local memory support, relative indexing into an array can compile to a chain of conditionals because there are no relative memory operations on arrays.

Some register-only tree traversal strategies developed the introduction of local and global memory write operations are appropriate for older hardware and may find renewed value for performance even on newer targets [Horn2007Stack] [Popov2007Stackless] .

Analytic and numerical roots can be mixed within a scene, as we did for the ground plane. Our simple ray tracer produces depth buffer values compatible with OpenGL's triangle rasterization, so implicit surfaces may be mixed freely with all primary ray strategies.

11. Some Online Examples in GLSL

Educational

I wrote these three heavily-documented examples to demonstrate how to set up an analytic and a ray marching renderer. These are on the Shadertoy website, which uses WebGL, a limited version of the full OpenGL supported by G3D. So, the code is a little less efficient and clear than the best case, but has the advantage of running in most current web browsers.

© 2013 Morgan McGuire
Analytic Tracer

© 2013 Morgan McGuire
Mandelbulb Explained

© 2013 Morgan McGuire
Legos. (Lego is a trademark of The Lego Group)

Aspirational

These three examples are by Iñigo Quilez, a graphics professional who's worked at Pixar and Oculus Story Studio. He's long contributed to the demoscene and online graphics community. These are three of the most impressive demos on the Shadertoy website and good (if a bit hard for the novice to parse) examples of how to use the implicit surface techniques described in this document.

© 2014 Iñigo Quilez
Dolphin.

© 2013 Iñigo Quilez
Mike.

© 2013 Iñigo Quilez
Elevated.

12. Further Reading

Implicit surfaces were first rendered using numerical root finding by Blinn [1982] . They've since been advanced significantly through research by many others. Gomes et al. give a good survey in their book [2009] .

The academic graphics community has largely worked in parallel with the demoscene community, which has been ray tracing these in real time with similar techniques since at least the early 2000's, but which has only recently begun to document and share their methods widely. Quilez has been a leading voice in spreading information about demoscene techniques through his personal website, the Shadertoy website, and recently some public presentations [Quilez20084k] and online videos [Quilez2014Video] . Swoboda's GDC 2012 presentation [2012] overviews the material from this document and then describes some of the performance issues on modern GPUs.

I've simplified some of the details of sphere tracing and distance estimators in this explanation. See Hart's SIGGRAPH course notes [1993] or Liktor's survey paper [2008] for more details on the mathematical constraints and ways of choosing optimal step sizes based on derivatives. See Keinert et al.'s [2014] paper for the current state of the art.

One of the most interesting applications of numerical methods for ray intersection is rendering fractals, which by definition have difficult surfaces to intersect (since they have unbounded area!) [Hart1989Fractal] .

An alternative approach for rendering surfaces that do not admit analytic intersection solutions is to convert them to an approximation that does, for example, a triangle mesh. A popular method for doing so is Lorensen's and Cline's marching cubes algorithm [1987] . Conversely, it is often desirable to convert an existing triangle mesh into a signed distance function for integration with a ray marcher (e.g., to achieve warping and smooth blending). Akleman and Chen [1999] give one implementation of such an operation. Swoboda [2012] discusses both triangle-to-implicit conversion and reverse in the context of modern GPUs.

The particular form of implicit surface defined by distance from a point set is very popular for modeling fluids. Fedkiw worked extensively with these and the book that he coauthored is a good reference [Osher2003LevelSet] . There's an entire subfield of rendering called point-based graphics that directly ray traces such representations, sometimes using implicit surface methods.

References

[Akleman1999Distance]

Generalized Distance Functions Ergun Akleman and Jianer Chen in 1999 Shape Modeling International, p. 72-79, March1, 1999. Official URL: http://doi.ieeecomputersociety.org/10.1109/SMA.1999.749326 10.1109/SMA.1999.749326

[Blinn1982Blobs]

A Generalization of Algebraic Surface Drawing James F. Blinn ACM Transactions on Graphics 1:3, p. 235-256, ACM, New York, NY, USA, July1982. Official URL: http://doi.acm.org/10.1145/357306.357310 10.1145/357306.357310

[Gomes2009Implicit]

Implicit Curves and Surfaces: Mathematics, Data Structures and Algorithms Abel Gomes, Irina Voiculescu, Joaquim Jorge, Brian Wyvill, and Callum Galbraith Springer Publishing Company, Incorporated, 2009.

[Hart1989Fractal]

Ray Tracing Deterministic 3-D Fractals J. C. Hart, D. J. Sandin, and L. H. Kauffman in Proceedings of the 16th Annual Conference on Computer Graphics and Interactive Techniques, p. 289-296, ACM, New York, NY, USA, 1989. Official URL: http://doi.acm.org/10.1145/74333.74363 Free URL: http://graphics.cs.illinois.edu/sites/default/files/rtqjs.pdf 10.1145/74333.74363

[Hart1993Course]

Ray Tracing Implicit Surfaces John C. Hart p. 1-16, SIGGRAPH'93 Course Notes: Design, Visualization and Animation of Implicit Surfaces, 1993. Free URL: http://graphics.cs.illinois.edu/sites/default/files/rtis-tr.pdf

[Hart1996Sphere]

Sphere Tracing: A Geometric Method for the Antialiased Ray Tracing of Implicit Surfaces John C. Hart The Visual Computer 12:10, p. 527-545, Springer, 1996. Free URL: http://bit.ly/1KJcnsO

[Hoffman1985Blend]

Automatic Surface Generation in Computer Aided Design C. Hoffman and J. Hopcroft The Visual Computer, p. 92-100, 1985. Free URL: http://www.cs.purdue.edu/homes/cmh/distribution/papers/Geometry/geo1.pdf

[Horn2007Stack]

Interactive K-D Tree GPU Raytracing Daniel Reiter Horn, Jeremy Sugerman, Mike Houston, and Pat Hanrahan in Proceedings of the 2007 Symposium on Interactive 3D Graphics and Games, p. 167-174, ACM, New York, NY, USA, 2007. Official URL: http://doi.acm.org/10.1145/1230100.1230129 10.1145/1230100.1230129

[Keinert2014Sphere]

Enhanced Sphere Tracing Benjamin Keinert, Henry Schäafer, Johann Korndörfer, Urs Ganse, and Marc Stamminger in Proceedings of Smart Tools & Apps for Graphics, Andrea Giachetti (ed.), p. 1-8, EuroGraphics, 2014. Official URL: http://diglib.eg.org/EG/DL/LocalChapterEvents/ItalChap/STAG2014/001-008.pdf Free URL: http://lgdv.cs.fau.de/get/2234

[Liktor2008Implicit]

Ray Tracing Implicit Surfaces on the GPU Gábor Liktor Computer Graphics and Geometry Journal 10:3, 2008. Free URL: http://www.cescg.org/CESCG-2008/papers/TUBudapest-Liktor-Gabor.pdf

[Lorensen1987MarchingCubes]

Marching Cubes: A High Resolution 3D Surface Construction Algorithm William E. Lorensen and Harvey E. Cline in Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques, p. 163-169, ACM, New York, NY, USA, 1987. Official URL: http://doi.acm.org/10.1145/37401.37422 Free URL: https://www.cs.virginia.edu/johntran/GLunch/marchingcubes.pdf 10.1145/37401.37422

[Osher2003LevelSet]

Level Set Methods and Dynamic Implicit Surfaces Stanley Osher and Ronald P. Fedkiw Springer, New York, N.Y., 2003. Official URL: http://opac.inria.fr/record=b1099358

[Popov2007Stackless]

Stackless KD-Tree Traversal for High Performance GPU Ray Tracing Stefan Popov, Johannes Gunther, Hans-Peter Seidel, and Philipp Slusallek Computer Graphics Forum 26:3, p. 415-424, 2007. Official URL: http://dblp.uni-trier.de/db/journals/cgf/cgf26.html#PopovGSS07

[Quilez20084k]

Rendering Worlds With Two Triangles With Raytracing on the GPU in 4096 Bytes Iñigo Quilez Talk at NVScene, August22, 2008. Free URL: http://www.iquilezles.org/www/material/nvscene2008/rwwtt.pdf

[Quilez2008Distance]

Modeling With Distance Functions Iñigo Quilez 2008. Official URL: http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm

[Quilez2014Video]

Formulanimations Tutorial :: The Principles of Painting With Maths Iñigo Quilez YouTube Video, July14, 2014. Free URL: https://www.youtube.com/watch?v=0ifChJ0nJfM

[Swoboda2012Procedural]

Advanced Procedural Rendering in DirectX 11 Matt Swoboda Talk at the Game Developers Conference, 2012. Free URL: http://directtovideo.wordpress.com/2012/03/15/get-my-slides-from-gdc2012/